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Abstract 

Liapunov  functions  are  constructed  for  nonlinear  systems  of  ordinary  differential  equations 
whose  linearized  system  at  an  equilibrium  point  possesses  either  a  simple  zero  eigenvalue  or 
a  complex  conjugate  pair  of  simple,  pure  imaginary  eigenvalues.  The  construction  is  explicit, 
and  yields  parametrized  families  of  Liapunov  functions  for  such  systems.  In  the  case  of  a  zero 
eigenvalue,  the  Liapunov  functions  contain  quadratic  and  cubic  terms  in  the  state.  Quartic 
terms  appear  as  well  for  the  case  of  a  pair  of  pure  imaginary  eigenvalues.  Predictions  of  local 
asymptotic  stability  using  these  Liapunov  functions  are  shown  to  coincide  with  those  of  pertinent 
bifurcation-theoretic  calculations.  The  development  of  the  paper  is  carried  out  using  elementary 
properties  of  multilinear  functions.  The  Liapunov  function  families  thus  obtained  are  amenable 
to  symbolic  computer  coding. 
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1  INTRODUCTION 


In  this  paper,  we  construct  families  of  Liapunov  functions  useful  in  assessing  the  asymptotic  stability 
of  critical  equilibrium  points  of  a  class  of  systems 

i  =  f(x)  (1) 

where  x  G  IRn  and  /  is  at  least  four  times  continuously  differentiable.  Throughout  the  presentation, 
we  let  the  origin  x  —  0  be  the  equilibrium  point  of  interest  of  this  system  of  ordinary  differential 
equations.  Critical  cases  in  the  study  of  stability  of  the  origin  of  Eq.  (1)  are  those  in  which  the 
Jacobian  matrix  D/( 0)  possesses  at  least  one  eigenvalue  with  zero  real  part,  and  no  eigenvalues 
with  positive  real  part.  In  these  situations,  it  is  not  possible  to  ascertain  whether  or  not  the  origin 
of  (1)  is  locally  asymptotically  stable  solely  from  the  linearized  system  x  =  ( Df(0))x . 

We  focus  on  two  specific  critical  cases  in  the  stability  analysis  of  Eq.  (1).  These  correspond 
to  the  Jacobian  matrix  of  (1)  at  the  origin  possessing  either  a  simple  zero  eigenvalue  or  a  complex 
conjugate  pair  of  simple,  pure  imaginary  eigenvalues.  The  associated  hypotheses  ((S)  and  (H) 
below,  respectively)  are  introduced  next.  The  techniques  of  this  paper  can  also  be  applied  to  other 
critical  cases,  such  as  those  involving  multiple  critical  eigenvalues  [7]. 

In  essence,  the  paper  relies  on  two  basic  tools  in  achieving  its  goals.  First,  the  notation  of 
multilinear  functions  is  adhered  to  throughout  the  paper,  in  denoting  terms  in  the  Taylor  series 
expansions  of  the  nonlinear  system  of  interest  as  well  as  Liapunov  function  candidates  and  their 
derivatives.  Second,  a  new  result  on  local  definiteness  of  a  class  of  scalar  bivariate  functions 
is  introduced.  This  appears  as  Lemma  1  in  Section  4.  Lemma  1  is  indeed  key  to  the  subsequent 
constructions  of  Liapunov  functions  in  Cases  (S)  and  (II).  These  constructions  are  rather  mechanical 
given  Lemma  1,  albeit  somewhat  tedious. 

The  first  critical  case  of  interest  in  this  work  is  characterized  by  the  occurrence  of  a  simple 
zero  eigenvalue  of  Df{ 0),  with  all  remaining  eigenvalues  having  strictly  negative  real  parts.  As 
discussed  in  [2],  stability  of  the  origin  in  this  situation  is  closely  related  to  the  stability  of  bifurcated 
equilibrium  points  in  smooth  parametrizations  of  Eq.  (1).  Because  of  this  connection  to  stationary 
(or  static)  bifurcation,  this  critical  case  will  be  referred  to  here  as  “Case  (S).” 

(S)  The  Jacobian  0/(0)  possesses  a  simple  zero  eigenvalue,  with  all  other  eigenvalues  in  the  open 
left  half  of  the  complex  plane. 
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In  the  second  critical  case  of  interest  here,  £)/( 0)  is  assumed  to  possess  a  single  complex  con¬ 
jugate  pair  of  simple,  pure  imaginary  eigenvalues,  with  the  remaining  eigenvalues  lying  in  the  open 
left  half  complex  plane.  Reference  [1]  discusses  the  relationship  of  stability  of  the  origin  in  this  case 
with  stability  of  bifurcated  periodic  solutions  of  smooth  parametrizations  of  (1).  Hopf  bifurcation 
to  periodic  solutions  occurs  for  generic  such  parametrizations  of  Eq.  (1)  under  these  circumstances. 
Because  of  this  connection  to  Hopf  bifurcation,  this  critical  case  will  be  referred  to  as  “Case  (H).” 

(H)  The  Jacobian  D/( 0)  possesses  a  complex  conjugate  pair  of  simple,  pure  imaginary  eigenvalues, 
with  all  other  eigenvalues  in  the  open  left  half  of  the  complex  plane. 

Recently,  there  has  been  significant  interest  in  feedback  stabilization  of  nonlinear  systems  in 
critical  cases  (see,  e.g.,  the  review  paper  [19]  and  references  therein).  This  has  yielded  various 
existence  and  synthesis  results  on  stabilizability  by  either  smooth  or  continuous  feedback.  The  main 
contribution  of  the  present  paper  to  this  body  of  work  is  the  construction  of  new  Liapunov  functions 
for  such  systems,  in  the  two  critical  cases  (S)  and  (II).  Liapunov  functions  facilitate  estimation  of 
the  domain  of  attraction  of  a  stable  equilibrium  point,  and  as  such  can  serve  to  quantify  the 
efficacy  of  a  given  control  design.  Liapunov  functions  can  be  used  to  define  performance  indices 
in  optimization-based  feedback  control  design  of  nonlinear  systems.  Such  performance  indices  can 
involve  estimates  of  the  achieved  domain  of  attraction  and  measures  of  the  adequacy  of  the  transient 
response.  The  Liapunov  functions  derived  here  for  Cases  (S)  and  (H)  are  given  explicitly  in  terms 
of  the  system’s  dynamics,  and  are  amenable  to  symbolic  computer  coding.  As  a  by-product  of  our 
results,  known  formulae  for  testing  stability  in  the  critical  cases  (S)  and  (H)  (so-called  bifurcation 
formulae )  are  found  to  follow  easily  from  the  Liapunov  functions  we  obtain.  Specifically,  in  Case 
(S)  it  is  known  that  the  vanishing  of  a  coefficient  /?i  and  the  negativity  of  a  coefficient  /?2  together 
imply  local  asymptotic  stability  of  the  origin.  Expressions  for  these  coefficients  appear  in  Eqs. 
(51)  and  (52)  below,  respectively.  Similarly,  in  Case  (H),  local  asymptotic  stability  of  the  origin  is 
implied  by  the  negativity  of  a  coefficient  /?2,  an  expression  for  which  appears  in  Eq.  (68)  below.  The 
bifurcation  formulae  (51),  (52),  (68)  are  expressions  based  on  the  linear,  quadratic  and  cubic  terms 
in  the  Taylor  expansion  of  f(x )  near  0.  The  stability  information  obtained  from  these  formulae  is 
local,  i.e.,  valid  in  a  sufficiently  small  neighborhood  of  the  origin.  Thus,  these  formulae  alone  do 
not  yield  analytical  performance  indices  of  the  type  just  alluded  to. 
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The  asymptotic  stability  of  nonlinear  systems  in  critical  cases  has  received  significant  attention 
in  the  literature  (e.g.,  [3],  [11],  [13]-[15] ,  [17],  [20]).  Although  some  of  these  works  have  employed 
Liapunov  stability  analysis,  the  Liapunov  functions  used  have  generally  been  defined  only  implicitly. 
In  some  cases,  this  is  linked  to  the  use  of  implicitly  defined  nonlinear  coordinate  transformations 
to  lower  dimensional  problems.  Implicitly  defined  Liapunov  functions  suffice  when  the  goal  of  the 
analysis  is  limited  to  deriving  sufficient  conditions  for  local  asymptotic  stability.  For  instance,  one 
result  of  Mees  and  Chua  [17]  gives  a  Liapunov  function  for  second  order  systems  (1)  satisfying 
(H).  Implications  of  this  result  for  higher  dimensional  models  (1)  follow  from  the  Center  Manifold 
Theorem  (cf.  [17], [5]).  In  the  present  paper,  we  give  an  explicit  construction  of  families  of  Liapunov 
functions  for  critical  nonlinear  systems  satisfying  either  hypothesis  (S)  or  (H)  which  apply  directly 
to  the  given  n-dimensional  system  description  (1). 

This  paper  is  organized  as  follows.  In  Section  2,  pertinent  results  on  multilinear  functions 
are  given.  The  set-up  for  construction  of  Liapunov  functions  for  systems  (1)  is  formulated  using 
multilinear  function  notation  in  Section  3.  A  lemma  giving  sufficient  conditions  for  local  definiteness 
of  a  class  of  scalar  bivariate  functions  is  presented  in  Section  4.  Section  5  contains  a  result  on 
solutions  of  Liapunov  matrix  equations  for  a  coefficient  matrix  with  a  zero  eigenvalue  or  a  pair  of 
pure  imaginary  eigenvalues.  The  main  results  of  the  paper  appear  in  Sections  6  and  7.  Section 
6  contains  an  explicit  construction  of  a  family  of  Liapunov  functions  for  Case  (S),  and  Section  7 
contains  an  analogous  construction  for  Case  (H).  Conclusions  are  collected  in  Section  8.  Appendix 
A  contains  a  discussion  of  linear  algebraic  equations  with  a  singular  coefficient  matrix. 

Notation.  In  what  follows,  IRn  denotes  the  space  of  n-dimensional  column  vectors  having  real 
entries,  while  47n  denotes  the  space  of  n-dimensional  column  vectors  with  complex  entries.  The 
complex  conjugate  of  a  quantity  (scalar,  vector,  or  matrix)  a  is  denoted  by  a.  The  transpose  of  a 
vector  or  matrix  a  is  denoted  aT .  For  a  vector  space  V,  denote  by  ( V)k  the  vector  space  obtained  as 
the  fc-tuple  product  V  x  ■  ■  ■  x  V .  The  Jacobian  derivative  of  a  function  <j>  is  denoted  D<\>.  The  norm 
of  a  vector  x  £  IRn  will  be  denoted  |x|,  and  the  same  notation  will  apply  to  any  compatible  matrix 
norm.  Denote  by  r  (resp.  /)  the  right  column  (resp.  left  row)  eigenvector  of  D/( 0)  corresponding 
to  the  critical  eigenvalue  0  (Case  (S))  or  iwc  (Case  (H)).  For  consistency  with  previous  literature 
[9],  [1],  [2],  the  first  component  of  r  is  set  to  unity,  and  /  is  then  chosen  subject  to  the  normalization 
lr  =  1.  (Ensuring  that  the  first  component  of  r  is  nonzero  in  some  cases  requires  a  reordering  of  the 
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elements  of  x.)  Finally,  given  a  matrix  A  with  a  simple  zero  eigenvalue,  A~  denotes  the  restricted 
inverse  of  A  given  by  A~  :=  (ATA  +  lTl)~lAT  (Eq.  (A. 5)  of  Appendix  A),  where  /  is  the  left 
eigenvector  of  A  corresponding  to  the  zero  eigenvalue,  normalized  in  the  manner  just  discussed. 

2  RESULTS  ON  MULTILINEAR  FUNCTIONS 

Multivariable  Taylor  series  can  be  conveniently  represented  in  terms  of  multilinear  functions.  We 
shall  employ  multilinear  functions  in  representing  Taylor  series  expansions  both  for  the  vector  field 
/(x)  of  Eq.  (1),  and  for  the  Liapunov  functions  whose  construction  is  the  main  purpose  of  this 
work.  In  this  section,  we  present  several  useful  facts  pertaining  to  multilinear  functions. 

2.1  Multilinear  Functions 

Multilinear  functions  may  be  defined  as  follows. 

Definition  1.  Let  Vj ,  V2,. . .  ,14  and  W  be  vector  spaces  over  the  same  field.  A  map  ip  :  Vi  x  V2  x 
. . .  x  14  1— *  W  is  said  to  be  multilinear  (or  k-linear )  if  it  is  linear  in  each  of  its  variables.  That  is 
[4,  p.  76],  for  arbitrary  vl,vl  £  Vj,  i  =  1, . . .  ,k,  and  for  arbitrary  scalars  a,  a,  we  have 

ip(vx , . . .  ,av'  +  av\...,vk)  =  aip^1 , . . . ,  v\ . . . ,  vk)  +  aip^1 , . . . ,  v* , . . . ,  vk).  (2) 


We  refer  to  k  as  the  degree  of  the  multilinear  function  ip.  In  particular,  multilinear  functions  of 
degree  two,  three  and  four  are  referred  to  as  bilinear ,  trilinear  and  tetralinear  functions,  respectively. 

We  shall  in  the  sequel  deal  exclusively  with  multilinear  functions  ip  whose  domain  is  the  product 
space  of  k  identical  vector  spaces  Vi  =  V2  =  ■  ■  •  =  Vk  =  V.  For  such  multilinear  functions,  we  have 
the  following  notion  of  symmetry. 

Definition  2.  A  fc-linear  function  ip  :  V  x  V  x  ...  x  V  w  IT  is  symmetric  if,  for  any  vl  £  V, 
i  —  1, . . . ,  k,  the  vector 

ip(v1,v2,...,vk)  (3) 

is  invariant  under  arbitrary  permutations  of  the  argument  vectors  vl .  ■ 

With  an  arbitrary  multilinear  function  ip,  we  associate  a  symmetric  multilinear  function  Sym  ip 
resulting  from  the  following  simple  device,  known  as  the  symmetrization  operation  [4,  pp.  88-89],  [6, 
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p  -  2  2] .  Given  a  multilinear  function  ^(x1,  x2, . . . ,  xk),  define  anew  (symmetric)  multilinear  function 
Sym  ip  as  follows: 

Symip(x1,x2,...,xk)  :=  ^  J2  H*'1  >  xh , . . . ,  xik),  (4) 

(l'l  ,t2,...,u) 

where  the  sum  is  taken  over  the  k\  permutations  of  the  integers  1,2,  ...,k. 

2.2  Jacobian  of  a  Polar  Form 

Recall  that  a  function  <j>  :  IRn  i— *  lRm  is  homogeneous  of  degree  k  ( k  an  integer)  if  (f>(ax)  —  ak<p(x) 
for  any  a  E  IR ,  x  E  IRn ■  Functions  which  are  polynomic  in  x  and  homogeneous  of  degree  k, 
k  >  1,  arise  naturally  in  Taylor  series  calculations.  This  section  pertains  to  such  functions,  and 
their  relationship  to  multilinear  functions. 

Definition  3.  Given  a  multilinear  function  ip  :  ( lRn)k  IRm ,  define  ipp,  the  polar  form  induced 
by  ip  [6,  p.  22],  as 

ipp(x)  ip(x,x,...,x).  (5) 


In  the  cases  k  =  2,3  and  4,  ipp  is  referred  to  as  a  quadratic,  cubic  and  quartic  form,  respectively. 
Evaluating  the  Jacobian  of  a  polar  form  ipp  is  facilitated  by  the  following  result,  which  is  a 
variant  of  Euler’s  Theorem  on  Homogeneous  Functions  (see,  e.g.,  [18,  p.  120]).  The  symmetrization 
operation  (4)  is  useful  in  exhibiting  a  symmetric  multilinear  function  ip(xl ,  x2, . . . ,  xk )  that  induces 
the  polar  form  ipp,  a  preliminary  step  in  applying  Proposition  1. 

Proposition  1.  (Jacobian  of  a  Polar  Form)  Let  ip  :  ( IRn)k  i— ►  IRm  be  a  symmetric  ^-linear 
function.  With  ipP{x)  :=  tp(x,  x, . . .  ,x),  and  v  an  arbitrary  vector  in  IRn , 

(Dipp(x))v  =  kip(x,x, . . .  ,x,v).  (6) 

Proof.  Fix  x  E  IRn ,  and  let  6x  E  IRn  be  arbitrary.  Then,  invoking  first  the  multilinearity  and 
then  the  symmetry  of  ip,  we  have 

fp (~T^)  =  f  (~  T  *  T  -  -  • ,  *  T  Sx) 

=  ip{x,x, . . .  ,*)  +  ip(6x,x,  x,...,x)+  ip(x,  8x,  x,  +  ip(x,x,  ...,x,6x) 
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(7) 


+0(|^|2) 

=  tpp(x)  +  k\j>(x,  X, . . . ,  X,  8x)  +  0(|<5i|2). 

On  the  other  hand,  by  Taylor’s  formula, 

t/)p(x  +  8x)  =  xl>p{x)  +  ( Dipp(x))8x  +  0{\8x\2).  (8) 

Comparing  Eqs.  (7)  and  (8),  we  obtain  the  desired  result.  ■ 

2.3  Coordinate  Representation  of  Scalar  Multilinear  Functions 

In  this  subsection,  we  state  a  useful  representation  result  for  scalar  multilinear  functions  ip  : 
(IRn)k  i — ?  IRm .  The  representation  rests  upon  a  choice  of  basis  (“coordinates”)  for  IRn .  Thus, 
let  {r1 ,  r2, .  . . ,  r”}  be  a  basis  for  IRn .  By  a  standard  result  [4,  Proposition  3-6.1],  to  this  basis  there 
corresponds  a  unique  dual  basis  which  we  may  view  as  consisting  of  row  vectors  l1  ,l2, . . .  ,ln  such 
that 

/V  =  Sii  (9) 

for  i,j  =  1, . . . ,  n.  Here,  8{j  is  the  Kronecker  delta  symbol: 

1  if  i  —  h 
0  if  i  ^  j 

for  i,j  —  1 , ,n.  Eq.  (9)  will  be  referred  to  as  the  biorthonormality  property  of  the  vectors  ll ,  rJ . 

Proposition  2  below  provides  a  convenient  representation  for  scalar  symmetric  multilinear  func¬ 
tions  on  ( IRn)k  in  terms  of  the  dual  basis  vectors  ll  and  a  set  of  “structural  coefficients.”  The  result 
will  be  applied  in  the  next  section,  yielding  a  representation  of  Liapunov  function  candidates. 
Proposition  2.  (Coordinate  Representation  of  Scalar  Multilinear  Functions)  Any  sym¬ 
metric  ^-linear  function  :  (Zfi!”)^  h— ►  IR  can  be  written  as 

=  .(rVxffe*2)...  (/•'»**)  (ii) 

where  the  (£-tuple)  sum  is  taken  over  all  i\ ,  *2, . . . ,  ik,  and  where  the  structural  coefficients 
are  symmetric  with  respect  to  all  permutations  of  *i,  *2,  •  •  • ,  ik- 
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Proof.  This  result  is  a  special  case  of  [16,  Thm.  1.2],  However,  we  sketch  a  rather  straightforward 
proof  for  the  sake  of  completeness.  By  [4,  Thm.  2.12.2],  a  scalar  multilinear  function  is  determined 
by  its  values  when  evaluated  at  all  combinations  of  basis  vectors  as  arguments.  The  formula  (11) 
for  ip  is  clearly  that  of  a  ^-linear  function.  Moreover,  it  follows  from  the  biorthonormality  property 
(9)  of  basis  vectors  rl  and  dual  basis  vectors  P  that,  by  appropriate  assignment  of  the  structural 
coefficients  ipi1i2...ik,  any  set  of  such  values  may  be  achieved.  Hence,  the  representation  above  is 
sufficiently  general  to  accomodate  any  scalar  multilinear  function  ip.  ■ 

2.4  Complexification  of  Real  Multilinear  Functions 

It  is  sometimes  convenient  to  evaluate  a  (real)  multilinear  function  ip  :  ( IRn)k  i— ►  IRm  for  argument 
vectors  in  47”.  This  is  done  simply  by  evaluating  the  value  of  ip  as  if  the  argument  vectors  were  in 
IRn,  using  a  representation  of  ip  such  as  Eq.  (11)  above.  This  process  is  the  complexification  of  ip. 

Our  use  of  the  complexification  device  is  relegated  to  Section  7,  in  the  construction  of  Liapunov 
functions  for  Eq.  (1)  under  hypothesis  (H).  The  following  observation  will  be  important  in  ensuring 
that  the  constructed  Liapunov  functions  are  indeed  real- valued. 

Proposition  3.  (Test  for  Realness  of  Multilinear  Functions)  Let  ip  denote  a  symmetric 
^-linear  function  ip  :  ( 47")fc  ^  <Sm .  The  image  of  ( IRn)k  under  the  map  ip  is  IRm  if  and  only  if 

ip(xl,x2, . . .  ,xk)  =  ip(xl,x2, . . .  ,xk )  (12) 

for  all  vectors  x1 ,x2, ...  ,xk  £  47”. 

Proof.  The  “ if  ”  part  is  automatic.  An  induction  proof  is  now  sketched  for  the  “ only  if”  part. 
Let  j  denote  the  number  of  argument  vectors  x‘  £  47”  that  are  not  also  in  IRn .  That  (12)  holds 
when  j  =  0  is  obvious.  Also,  if  (12)  holds  for  some  j  =  jo  <  then  it  is  a  simple  exercise  to  verify 
that  it  also  holds  for  j  =  j0  -f  1.  ■ 

3  REPRESENTATION  OF  LIAPUNOV  FUNCTION 
CANDIDATES 

The  Liapunov  functions  constructed  in  this  paper  consist  of  quadratic,  cubic,  and,  in  Case  (II), 
quartic  terms  in  the  state.  The  fact  that  quadratic  Liapunov  functions  exist  in  some  cases  might 
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lead  one  to  ask  whether  or  not  there  is  a  need  for  inclusion  of  terms  of  order  three  or  higher. 
Therefore,  we  address  this  issue  at  the  outset  by  giving  two  examples,  one  in  Case  (S)  and  one  in 
Case  (H),  in  which  no  quadratic  Liapunov  function  exists,  although  the  origin  is  asymptotically 
stable.  Liapunov  functions  proving  stability  for  these  two  examples  follow  easily  from  the  algorithms 
given  in  Sections  6  and  7. 

Example  1.  The  system 

ii  =  x\x2  -f  xf 
x2  =  — *2  ~  2x2 


satisfies  hypothesis  (S),  and  the  origin  is  asymptotically  stable.  However,  it  is  straightforward 
to  check  that  no  quadratic  Liapunov  function  exists  for  the  origin.  Liapunov  functions  can  be 
constructed  containing  quadratic  and  cubic  terms.  One  such  Liapunov  function  is 

x2  7 

V(xi,  x2)  =  x2  +  +  -x2x2. 


Example  2.  The  system 

xi  =  -x2  +  x2 

x2  —  x\  T  X2 

satisfies  hypothesis  (H).  Moreover,  the  origin  is  locally  asymptotically  stable,  as  can  be  verified  by 
computing  the  coefficient  /?2  using  Eqs.  (66)-(68)  below  and  checking  that  /?2  <  0.  However,  it 
is  straightforward  to  check  that  no  quadratic  Liapunov  function  exists  for  the  origin.  Liapunov 
functions  can  be  constructed  containing  quadratic,  cubic  and  quartic  terms,  using  the  algorithm 
of  Section  7  or,  since  this  is  a  second-order  system,  the  results  of  Mees  and  Chua  [17].  One  such 
Liapunov  function,  obtained  using  the  results  of  [17],  is 

V(x!,x2)  =  ^(xl  +  x22)  -  xIx2  +  y--Y 

XlX2  XjX2  x\ 

+~ 4  4  2~ 
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Indeed,  the  time  derivative  of  V  along  trajectories  of  this  system  is  —  |(x2  +  x2)2 .  ■ 

Eq.  (1)  may  be  rewritten,  upon  Taylor  series  expansion  of  /(x),  in  the  form 

x  =  f  0) 

=  Ax  +  Q(x,x)  +  C(x,x,x)  H -  (13) 

Here,  A  :=  D/( 0)  and  Q(x,x),  C(x,x,x)  are  vector-valued  quadratic  and  cubic  forms,  with  the 
dots  denoting  higher  order  terms.  Without  loss  of  generality,  assume  that  Q(x,x)  is  induced  by 
a  symmetric  bilinear  function  Q(x1,x2)  and,  similarly,  that  C(x,x,x )  is  induced  by  a  symmetric 
trilinear  function  C(x* ,  x2,  x3). 

In  Case  (S)  (one  zero  eigenvalue),  we  shall  in  the  sequel  seek  Liapunov  functions  V(x)  consisting 
of  the  sum  of  a  quadratic  part  and  a  cubic  part,  viz. 

V(x)  =  xTVx  +  Q(x,x,x).  (14) 

It  is  natural  to  require  V  to  be  symmetric  and  positive  definite.  Similarly,  the  cubic  form  <J(x,  x,  x) 
is  induced  by  a  symmetric  trilinear  function  £/(x* ,  x2,  x3).  Since  the  term  Q(x,x,x)  is  dominated 
by  the  quadratic  term  xLVx ,  any  such  V(x)  will  indeed  be  locally  positive  definite. 

In  our  study  of  Case  (H)  (two  purely  imaginary  eigenvalues),  we  will  include  a  quartic  term 
T(x,x,x,x)  in  the  candidate  Liapunov  function  V(x)  in  addition  to  the  quadratic  and  cubic  terms 
presented  in  (2),  viz. 

V(x)  =  xTVx  +  Q(x,x,x)  +  T(x,x,x,x).  (15) 

We  of  course  ask  that  the  quartic  form  T (x,  x,x,x)  be  induced  by  a  symmetric  tetralinear  function 
T(xa ,  x2,  x3,  x4).  Certainly,  the  local  positive  definiteness  of  V(x)  with  V  >  0  remains  preserved 
under  inclusion  of  T (x,  x,  x,  x)  or  terms  of  still  higher  order. 

Next,  we  invoke  Proposition  2  for  the  cases  k  —  2,3  and  4,  obtaining  coordinate  representations 
of  the  bilinear,  trilinear,  and  tetralinear  functions  x1?Vx2,  t/(xx ,  x2,  x3),  and  T(x1 ,  x2,  x3,  x4), 
respectively. 

Consider  first  the  case  k  =  2.  Since  any  real  quadratic  form  x1  Vx2  is  determined  by  a  real 
symmetric  matrix  V ,  we  can  apply  Proposition  2  to  conclude  that  all  such  matrices  V  have  the 
form 

n  n 

*,=x:£viTp.  (is) 

«— i  j= i 
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where  7 =  7Tj,-  are  real  coeffiecients. 

The  following  representations  for  trilinear  functions  Q  and  tetralinear  functions  T  also  follow 
from  Proposition  2: 

so*1,*2,*3)  =  EEE^^1)^2)^*3),  (i7) 

i=  1  j= 1  fc=l 

T(^ ,  x2,  x3,  z4)  =  E  E  E  E  ^P(/^1)(Px2)(/fcx3)(/^4),  (18) 

*=1  j=l  A;=l  p=l 

respectively.  These  multilinear  functions  are  rendered  symmetric  by  imposing  the  condition  that 
the  values  of  the  structural  coefficients  (7 ,•_,-*)  and  ( r,j jtp)  do  not  depend  on  the  order  of  the  indices. 
The  simple  representations  above  for  the  bilinear,  trilinear  and  tetralinear  functions  appearing  in 
the  Liapunov  function  candidates  V  imply  that  the  construction  ofV  is  tantamount  to  specification 
of  the  structural  coefficients  irij ,  7^ ,  rtjkp . 

In  the  calculations  to  follow,  the  goal  is  to  obtain  sets  of  structural  coefficients  Ttij,Jijk,Tijkp 
which  ensure  the  local  negative  definiteness  of  V,  the  time  derivative  of  the  Liapunov  function 
candidate,  along  trajectories  of  (1).  Of  course,  this  will  only  be  possible  under  assumptions  on 
Eq.  (1)  which  guarantee  local  asymptotic  stability  of  the  origin.  In  the  next  section,  a  result  is 
presented  on  local  definiteness  of  a  class  of  bivariate  functions.  The  sufficient  conditions  for  local 
definiteness  provided  by  this  result  will  facilitate  a  systematic  derivation  of  local  Liapunov  functions 
and  conditions  for  local  asymptotic  stability. 

4  LOCAL  DEFINITENESS  OF  A  CLASS  OF  BIVARIATE 
FUNCTIONS 

We  now  introduce  an  interesting  lemma  which  will  prove  to  be  an  important  tool  in  exhibiting 
conditions  for  local  negative  definiteness  of  the  time  derivative  of  a  Liapunov  function  candidate. 
Lemma  1.  (Local  Definiteness  of  a  Class  of  Bivariate  Functions)  The  scalar  bivariate 
function 

6{u,  v)  =  a2ou2  +  am?;4 

OO  q  0  00  q  a  tc 

+auuv  +  a2iu  v  +  a30u  +  ai3tm  +  a22u  v  +  a3iu  v  +  a40u  +  0(|(u,  u)|  )  (19) 
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in  the  real  variables  u  and  v  is  locally  negative  definite  near  (u,v)  =  (0,0)  if  a20  <  0,  a0 4  <  0  and 
|c*i2 1  <  2^020004-  Here,  0(|(u,  i>)|5)  denotes  terms  of  fifth  and  higher  order  in  |(it,u)|. 

Remark  1.  (a)  In  applying  this  lemma,  for  simplicity  we  shall  only  use  a  special  form  of  the 
sufficient  condition  it  provides.  Namely,  S(u,v )  is  locally  negative  definite  if  a20  <  0,  a 04  <  0  and 
al2  =  0.  (b)  It  is  convenient  to  view  this  sufficient  condition  (i.e.,  that  of  part  (a))  as  follows. 

Consider  6(u,  v ),  a  scalar  polynomial  function  of  the  scalars  u  and  v,  for  which  the  leading  term  in 
6(u,  0)  is  a2 ou2,  and  the  leading  term  in  6(0, «)  is  004V4.  Then  the  assertion  is  that  6(u,  v)  is  locally 
negative  definite  when  two  basic  conditions  are  fulfilled:  First,  the  absence  of  the  two  terms  ao3f3, 
(112UV2]  Second,  the  local  negative  definiteness  of  the  univariate  functions  6(u,0)  and  6(0,  v)  (i.e., 
a 20  <  0  and  ao4  <  0,  respectively).  ■ 

Proof.  In  proving  Lemma  1,  we  neglect  the  terms  0(|(u,u)|5)  in  6(u,v),  since,  being  higher  order 
terms,  they  can  easily  be  incorporated  with  only  slight  modifications  in  the  analysis.  With  this 
understanding,  rewrite  6(ii,u)  in  the  form  of  a  quadratic  polynomial  in  u: 

S(u,  v )  =  (a2o  +  a2iu  +  a30u  +  a22«2  +  0311m  +  a40u2)u2  +  (a^v2  -f  a13u3)u  +  a04V4 

=  :  p(u,v)u2  +  (ai2v2  +  ai3v3)u  +  a04V4.  (20) 


Here, 

p(u,  v)  020  +  a2\v  +  a30u  +  a22u2  +  031111;  +  a.4002.  (21) 

Since  p(0,  0)  =  a20  <  0,  it  is  clear  that  there  is  an  ei  >  0  such  that  p(u,  v)  <  0  for 


M,  |u|  <  €4. 


(22) 


(One  could  easily  write  a  formula  for  such  an  ej.) 

The  leading  coefficient  p(u,v)  in  the  expression  (20)  for  6(u,u)  is  therefore  strictly  negative  for 
M,  M  <  fi-  Next,  rewrite  6(it,v)  as 


6(u,  u)  =  p{u,  u)[u  + 


(q12u2  +  Qi3i;3) 
2 p(u,v) 


]2  +  ?(«,«), 


(23) 


where 


q(u,v)  := 


4 p(u,  v ) 


{4a2o«04  —  al2 


+  4a3Qa04U  +  (4a2i  a04  —  2ai2ai3)v 
+  4ao404ou2  -)-  4ao4a3juv  +  (4oo4022  —  o23)u2}. 


(24) 
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Since  020  <  0,  004  <  0  and  |ai2 1  <  ^\/ci2Ga04 1  the  constant  term  in  the  expression  in  braces  in  Eq. 
(24),  namely  4020^04  —  ai2i  is  strictly  positive.  Hence,  there  is  an  >  0  such  that  the  expression 
in  braces  in  Eq.  (24)  is  strictly  positive  for 

M,  M  <  e2-  (25) 

Recalling  that  p(u,v )  <  0  for  |u|,  |w|  <  ex,  we  have  that  for  |«|,  |u|  <  e  :=  min(ex,e2),  g(u,w)  <  0, 
with  q{u,  w)  =  0  only  for  v  =  0  (see  Eq.  (24)).  Now  consider  the  implications  of  these  observations 
for  the  expression  (23)  for  <5(u,  v).  Clearly,  for  |u|,  |u|  <  e  and  v  ^  0,  the  fact  that  q(it,  v )  is  strictly 
negative  ensures  that  6[u,v)  <  0.  If,  on  the  other  hand,  |u|,  |n|  <  e  and  v  —  0,  then  8(u,v)  reduces 
to 

6(u,0)  =  p(u,0)u2 

<  0  (26) 

for  u  0.  Thus,  6(u,v)  is  indeed  locally  negative  definite  near  (0,0).  ■ 

5  CALCULATIONS  INVOLVING  THE  STABLE  SUBSPACE 

In  this  section,  we  define  the  stable  subspace  of  IRn  corresponding  to  the  Jacobian  matrix  A,  recall 
an  associated  orthogonality  property  from  [10],  and  employ  the  stable  subspace  concept  in  the 
choice  of  the  quadratic  term  xTVx  in  the  Liapunov  function  candidate  V(x)  (cf.  Eqs.  (14),  (15)). 
The  development  proceeds  for  Cases  (S)  and  (H)  in  parallel. 

Definition  4.  The  stable  subspace  of  IRn,  denoted  by  Es ,  is  the  span  of  the  eigenvectors  (and 
generalized  eigenvectors,  if  any)  of  A  corresponding  to  the  stable  eigenvalues  of  A.  ■ 

In  Case  (S),  any  vector  x  €  IRn  has  a  unique  representation  x  =  ar  +  w  where  a  is  a  real  scalar, 
r  is  the  right  eigenvector  of  A  corresponding  to  the  eigenvalue  0,  and  w  £  Es .  In  Case  (II),  any 
vector  x  £  IRn  has  a  unique  representation  x  =  ar  +  ar  +  w  where  a  is  a  complex  scalar,  r  is  the 
right  eigenvector  of  A  corresponding  to  the  eigenvalue  iu)c,  and  w  £  Es . 

The  following  property  is  well  known  (see,  e.g.,  [10,  Appendix  4.1]). 

Proposition  4.  (Orthogonality  of  Left  and  Right  Eigenvectors)  Let  la  and  denote  left 
and  right  eigenvectors,  respectively,  corresponding  to  eigenvalues  Xa  and  Xp  of  a  matrix  A.  Bather, 
or  both,  of  P  and  may  be  generalized  eigenvectors.  If  Xa  /  Xp,  then  Pr &  =  0.  Moreover, 
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the  subspace  of  all  column  vectors  nullified  by  la  is  precisely  the  span  of  all  right  eigenvectors  and 
generalized  right  eigenvectors  of  A  associated  with  eigenvalues  other  than  \a. 

Remark  2.  Proposition  4  implies  the  following  facts,  which  will  prove  useful  in  the  sequel.  As 
above,  let  l  denote  a  left  eigenvector  corresponding  to  the  critical  eigenvalue  0  (in  Case  (S))  or  iujc 
(in  Case  (H)).  Then  Iw  =  0  if  and  only  if  w  £  Es .  Moreover,  pw  =  0  for  a  row  vector  p  if  and 
only  if  p  £  span(/)  (Case  (S)),  or  p  £  span(9?e  /,$rm  l)  (Case  (II)).  Finally,  in  Case  (H),  we  have  If 
=  Jr  =  0.  ■ 

Since  the  Jacobian  matrix  A  —  DfJ 0)  has  part  of  its  spectrum  on  the  imaginary  axis,  it  is  not 
possible  to  choose  a  positive  definite  V  for  which  AT7>  -\-V  A  is  negative  definite.  However,  one  can 
ensure  that  the  latter  matrix  is  negative  definite  on  a  subspace  of  lRn ,  while  being  only  negative 
semidefinite  on  all  of  IRn .  A  method  for  achieving  this  is  given  next. 

Recall  that  r  denotes  the  eigenvector  of  A  corresponding  to  the  critical  eigenvalue  (0  in  Case 
(S),  ioJc  in  Case  (II)).  Note  that  r  £  (Vn  in  Case  H,  and  that  r  6  IRn  in  Case  (S).  Denote  by  Es  the 
subspace  of  IRn  spanned  by  the  eigenvectors  (and  generalized  eigenvectors,  if  any)  corresponding 
to  the  stable  eigenvalues  of  A  (in  either  Case  (S)  or  Case  (H)).  We  refer  to  Es  as  the  stable  subspace 
of  IRn. 

The  following  proposition  is  useful  in  selecting  the  quadratic  term  xTVx  in  the  Liapunov  func¬ 
tion  V  under  either  hypothesis  (S)  or  (H). 

Proposition  5.  (Liapunov  Matrix  Equation  on  Stable  Subspace)  Using  the  notation  above, 
and  under  either  hypothesis  (S)  or  (II),  there  exists  a  family  of  real  symmetric  n  x  n  matrices  n 
for  which 

(i)  nr  =  0,  (ii)  trTnu;  >  0,  and  (iii)  wT {ArW  +  n A)w  <  0  (27) 

for  all  w  £  Es ,  w  0- 

Proof.  The  proof  may  be  carried  out  in  two  steps.  In  the  first  step,  we  exhibit  a  choice  of 
coordinates  for  the  state  space  IRn  for  which  the  existence  of  matrices  n  is  transparent.  In  the 
second  step,  we  verify  that  the  existence  of  a  matrix  n  satisfying  (i)-(iii)  in  one  coordinate  system 
implies  the  existence  of  such  a  matrix  for  any  choice  of  coordinates.  Step  1.  Suppose,  then,  that 
the  state  x  of  system  (1)  is  expressed  with  respect  to  a  coordinate  basis  {r1,  r2, . . . ,  rn }  defined  as 
follows.  In  Case  (S),  take  r1  r  and  choose  the  remaining  basis  vectors  r‘  £  IRn ,  i  —  2 ,  ...,n 
such  that  spanjr2, . .  . ,  r")  =  Es .  Analogously,  in  Case  (II),  choose  r1  :=  3?e(r),  r2  :=  7rm(r),  and 
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let  r3, . . . ,  rn  satisfy  span{r3, . . . ,  rn}  =  Es .  For  such  a  coordinate  basis,  A  has  the  block  diagonal 
representation 


A=r  °). 

VO  As) 

Here,  As  is  a  real  square  stable  matrix  whose  eigenvalues  coincide  with  the  stable  eigenvalues  of  A, 
and  6  is  given  by 


6  = 


(29) 


in  Case  (S), 
in  Case  (H). 

(Note  that  9  is  a  scalar  for  Case  (S),  and  is  a  2  x  2  matrix  for  Case  (H).)  It  is  now  straightforward 
to  exhibit  a  matrix  n  satisfying  (i)-(iii).  Consider  the  matrix 


f  0 

\( 0 

V-^c 

0  J 

(°  °), 

Vo  n22; 


(30) 


where  n22  is  a  real  symmetric  positive  definite  matrix  of  dimension  (n  —  1)  in  Case  (S),  and 
dimension  (n  —  2)  in  Case  (H),  for  which  the  matrix 


A5Tn22  +  n22As 


(31) 


is  negative  definite.  The  existence  of  such  a  matrix  n22  is  clear,  since  As  is  stable.  Note  that,  for 
the  present  choice  of  coordinate  basis,  we  have  that  in  Case  (S)  the  first  component  of  any  vector 
w  £  Es  is  0,  and,  in  Case  (H),  the  first  two  components  of  w  are  0.  Also,  the  right  eigenvector  r 
is  given,  in  Case  (S)  and  Case  (H),  by  r  =  (1, 0, . . . ,  0)r  and  r  =  (1,  i,  0, . . . ,  0)T,  respectively.  The 
matrix  n  of  Eq.  (30)  is  now  easily  verified  to  satisfy  conditions  (i)-(iii).  Step  2.  Next,  we  show 
that  existence  of  a  matrix  n  satisfying  (i)-(iii)  in  one  set  of  coordinates  implies  existence  of  such 
a  matrix  for  any  set  of  coordinates.  Let  the  coordinate  change  be  determined  by  a  nonsingular 
transformation  matrix  E,  and  indicate  vectors  in  the  transformed  coordinates  by  a  hat  ('"): 


r  —  Er,  w  =  Eu>.  (32) 

In  the  new  coordinates,  A  has  the  representation 

A  =  EAXT1.  (33) 

It  is  straightforward  to  verify  that  the  matrix 

fi  =  (ET)-1ns-1  (34) 
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satisfies  conditions  (i)-(iii),  where  all  quantities  are  taken  in  the  new  coordinates.  (Note  that  the 
relationship  between  II  and  II  specified  in  Eq.  (34)  is  not  that  of  a  similarity  transformation.)  ■ 

As  a  corollary  to  Proposition  5,  we  have  the  following  result. 

Corollary  1.  Let  II  be  a  matrix  satisfying  (i)-(iii)  of  Proposition  5.  Then,  there  is  an  o  >  0  such 
that,  for  each  nonzero  w  £  Es , 

tnT(ATII  +  UA)w  <  —  a|u;|2.  (35) 

Proof.  This  follows  easily  using  the  proof  of  Proposition  5  and  the  standard  fact  that,  for  a 
negative  definite  real  matrix  Q ,  there  is  an  a  >  0  such  that  xTQx  <  —  a\x\2  for  all  x  £  lRn .  In  the 
present  setting,  the  role  of  Q  is  played  by  ASTU22  +  1122^*.  ■ 

In  our  construction  of  Liapunov  functions  for  Case  (S),  we  shall  employ  matrices  V  of  the  form 

V  -  n  +  lTl  (36) 

where  II  is  any  real  n  x  n  matrix  satisfying  the  conditions  of  Proposition  5.  For  Case  (II),  the 
matrices  V  will  be  taken  to  be  of  the  form 

V  =  H  +  lTI+lTl.  (37) 

It  is  not  difficult  to  check  that,  both  in  Case  (S)  and  Case  (H),  the  matrix  V  is  positive  definite, 
and  the  matrix  ATV  +VA  is  negative  se/m'definite,  but  negative  definite  when  restricted  to  the 
stable  subspace  Es .  This  class  of  positive  definite  matrices  V  will  be  sufficiently  general  for  the 
construction  of  Liapunov  function  candidates  in  Cases  (S)  and  (H). 

6  CONSTRUCTION  OF  LIAPUNOV  FUNCTIONS  IN  THE 
CASE  OF  ONE  ZERO  EIGENVALUE 

In  this  section,  we  construct  a  family  of  Liapunov  function  candidates  for  system  (1)  (equivalently, 
(13))  under  hypothesis  (S).  The  main  task  is  to  specify  the  matrix  V  and  the  cubic  form  Q(x,  x,  x ) 
appearing  in  the  expression  (14)  for  V(x).  In  the  foregoing,  we  have  constrained  the  matrix  V  to 
take  the  form  V  =  lT l  +  II  in  Case  (S),  where  II  is  any  real  symmetric  n  x  n  matrix  satisfying  the 
conditions  of  Proposition  5. 
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6.1  Conditions  for  V  <  0 


Using  Proposition  1  (Jacobian  of  Polar  Forms),  the  time  derivative  of  V(:r)  evaluated  along  trajec¬ 
tories  of  Eq.  (13)  can  be  written  as 

V(x)  =  xT(ATV  +VA)x 

+  2 QT(x,  x)Vx  4-  3G(x,  x,  Ax) 

+  2 CT(x,x,x)Vx  +  3 G(x,x,  Q(x,  *))+•••  (38) 

Recall  from  Section  5  the  notation  Es  for  the  stable  subspace  of  IRn,  i.e.,  the  (n—  l)-dimensional 
subspace  spanned  by  the  eigenvectors  (and  generalized  eigenvectors,  if  any)  corresponding  to  the 
stable  eigenvalues  of  A.  Using  the  representation  x  =  ar  +  w  (w  £  Es )  in  Eq.  (38);  recalling  that 
V  has  been  chosen,  by  Eq.  (36),  such  that  V  —  f  l  +  II  with  Hr  =  0;  invoking  the  fact  that  a 
multilinear  function  is  linear  in  each  argument;  and  collecting  terms  on  the  right  side  of  (38)  of  like 
order  in  |(a,  tn)|,  we  obtain  a  series  expansion 

VW:=[V](2)  +  [V](3)  +  [V](4)  +  .--  (39) 

where  the  integer  subscripts  denote  the  degree  of  the  corresponding  term  in  |(a,u))|,  and  the  dots 
denote  terms  of  fifth  and  higher  order  in  |(a,  tn)|.  Specifically,  the  terms  appearing  on  the  right  side 
of  Eq.  (39)  are  given  by 

[V](2)  =  wT(ATH  +  UA)w,  (40) 

[V](3)  =  2 a3lQ(r,  r ) 

+  3a2£/(r,  r,  Aw)  +  4a2lQ(r,  w)  +  2 a2QT  [r,  r)IIu) 

+  2 aIQ(w,  w)  +  4 aQT(r,  u;)nu;  4-  6at/(r,  w,  Aw) 

4-  2QT(w,  u;)IIu;  +  3(/(u>,  in,  Aw),  (41) 

[V](4)  =  a4{2lC(r,r,r)  +  3Q(r,r,Q(r,r))} 

+  2a3{3/C(r,  r,  w)  +  3 <7(r,  w,  Q(r,  r))  4-  3f7(r,  r,  Q(r,  u;))} 

+  3a2{2/C(r,  w,  w)  +  <y(r,  r,  Q(w,  w))  +  4Q{r,  u>,  Q(r,  w))  +  Q{w ,  w,  Q(r,  r))} 

4-  2 a{IC(w,  w ,  w)  -(-  3Q(r,  w,  Q(w,  iv))  +  3f7(u;,  w,  Q{r ,  w))} 

4-  2CT(w,  w,  u>)IIu;  +  3(7(tn,  w,  Q(w,  w )).  (42) 
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Note  that  condition  (iii)  on  IT  (cf.  Proposition  5)  implies  that  [V]j2)  <  0  for  w  6  E\w  ^  0. 
This  does  not  of  course  imply  that  V  is  locally  negative  definite,  only  that  it  is  locally  negative 
definite  on  the  subspace  Es . 

Lemma  1,  along  with  the  foregoing  computation  of  V,  allow  us  to  obtain  the  following  pre¬ 
liminary  statement  concerning  the  local  asymptotic  stability  of  the  origin  of  Eq.  (1).  Note  that 
condition  (SI)  of  the  next  proposition  is  a  known  necessary  condition  for  stability  for  systems  (1) 
possessing  a  simple  zero  eigenvalue  (see  for  instance  [2],  [10]). 

Proposition  6.  Under  hypothesis  (S),  the  origin  of  Eq.  (1)  is  locally  asymptotically  stable  if 
there  are  a  real  symmetric  nxn  matrix  II  satisfying  (i)-(iii)  of  Proposition  5,  and  a  symmetric  real 
trilinear  function  Q^x1 ,  x2 ,  x3),  for  which  the  following  three  conditions  hold: 

(51)  lQ(r,r)  =  0, 

(52)  3 Q(r,  r,  Aw)  +  4 lQ(r,  w)  +  2 Q1  (r,  r)IIu;  =  0  for  all  w  £  Es ,  and 

(53)  2 IC(r,  r,  r )  +  3 G(r,  r,  Q(r,  r))  <  0. 

Proof.  Let  conditions  (S1)-(S3)  of  the  Proposition  hold.  It  is  straightforward  to  write  an  upper 
bound  for  V(ar  +  w)  in  the  form  of  a  scalar  bivariate  Taylor  series.  To  facilitate  application  of 
Lemma  1,  we  employ  notation  consistent  with  Eq.  (19)  of  Lemma  1,  and  define  variables  u  :=  |u;| 
and  v  :=  |a|.  The  specialization  of  Lemma  1  given  in  Remark  1(a)  will  be  employed.  The  proof 
proceeds  in  two  steps.  The  first  step  consists  of  verifying  that  a03  =  0  and  a12  =  0  (see  Remark  1(a) 
following  Lemma  1).  In  the  second  step,  we  ascertain  that  a2 o  <  0  and  a0 4  <  0.  Step  1.  From  Eq. 
(41),  it  is  clear  that  the  n3-term  in  V  is  2/Q(r,  r)u3,  and  this  vanishes  by  virtue  of  condition  (SI). 
Thus,  the  upper  bound  for  V  will  naturally  be  absent  of  a  term  C103V3.  Using  Eq.  (41),  it  is  apparent 
that  the  uu2-term  in  the  upper  bound  will  vanish  if  3 Q(r,  r,  Aw)  +  4 lQ(r,  w )  +  2 QT(r,  r)IIu;  =  0  for 
each  w  £  Es .  This  latter  condition  is  precisely  (S2).  Step  2.  Eq.  (40)  and  Corollary  1  imply  that 
the  quadratic  terms  in  V  are  bounded  above  by  a  function  a2o«2,  where  a20  <  0.  Also,  Eq.  (42) 
(which  gives  the  quartic  terms  in  V)  and  assumption  (S3)  together  imply  that  004  <0.  ■ 

6.2  Algorithm  for  Construction  of  V  in  Case  (S) 

For  Proposition  6  to  be  useful  in  the  explicit  construction  of  Liapunov  functions  for  Eq.  (1),  a 
method  is  needed  facilitating  the  choice  of  a  matrix  II  and  a  trilinear  function  Q  for  which,  under 
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an  auxiliary  condition  guaranteeing  local  stability,  (S2)  and  (S3)  are  satisfied.  As  it  turns  out,  one 
can  first  choose  any  matrix  II  for  which  conditions  (i)-(iii)  of  Proposition  5  hold,  and  then  proceed 
to  construct  compatible  trilinear  functions  Q  satisfying  (S2),  (S3).  We  now  proceed  to  construct  a 
family  of  such  trilinear  functions,  using  the  representation  (17)  in  terms  of  the  associated  structural 
coefficients. 

The  general  representation  (17)  for  trilinear  functions  (/(x1,  x2,  x3)  assumes  a  specific  choice  of 
basis  for  IRn .  Let  {r1 ,  r2,  . . . ,  r”}  be  a  basis  for  IRn ,  obtained  by  setting  r1  :=  r  and  requiring 
that  r*  £  Es  for  i  =  2, . . .  ,  n.  Let  the  associated  dual  basis  (discussed  in  Section  2.3)  be  given 
by  {I1 ,  l2, . . . ,  /”}.  Recall  that  the  biorthonormality  property  (9)  holds,  i.e.,  that  Pr J  =  8{j,  the 
Kronecker  delta.  Moreover,  by  Proposition  4  and  the  fact  that  the  dual  basis  is  unique,  we  have 
that  I1  =  /,  the  left  eigenvector  of  A  associated  with  the  eigenvalue  zero  (recall  the  normalization 
Ir  —  1). 

By  (17),  the  trilinear  function  Q  can  be  represented  as  follows: 

g(x\x2,x3)  =  EEETiifc^*1)^*3)^*3).  (43) 

*=1  j—  1  k— 1 

Next,  this  representation  is  used  to  determine  all  the  trilinear  functions  Q  satisfying  condition 
(S2)  of  Proposition  6.  Using  (43)  and  the  biorthonormality  property  (9),  we  find  that  the  term 
Q(r,r,  Aw)  appearing  in  (S2)  is  given,  for  any  w  £  Es ,  by 

n  n  n 

Q(r,  r,  Aw)  = 

2—1 j= 1 k= 1 

n 

=  5>m(fMu;).  (44) 

i= 2 

From  Appendix  A,  the  linear  operator  A  is  invertible  when  restricted  to  the  subspace  Es .  The 
inverse  of  the  restricted  operator  is  A~  :=  (AT A  +  lTl)~1AT  (see  Appendix  A).  Thus,  knowledge 
of  a  vector  w  £  Es  is  tantamount  to  knowledge  of  w  :=  Aw,  also  a  vector  in  Es .  Moreover,  w  is 
given  in  terms  of  w  by  the  formula 

w  =  A~w 

=  ( ATA  +  lTl)-1ATw .  (45) 

Using  this  observation  and  the  representation  (44)  for  Q(r,r,  Aw),  we  have  that  condition  (S2) 
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(46) 


is  equivalent  to  the  statement  that,  for  all  w  £  E8 , 

n 

SY^injl'w)  +  4  lQ(r,  A~w)  +  2  QT(r,  r)UA~w  =  0. 

«= 2 

Since  Eq.  (46)  is  linear  in  w,  and  since  the  vectors  r\  i  =  2, . . . ,  n  form  a  basis  for  Es ,  condition 

(52)  is  tantamount  to  requiring  that  (46)  holds  for  w  —  r\  i  =  2, ...  ,n.  Substitution  of  these 
values  for  w  in  (46)  yields,  using  the  biorthonormality  property,  that,  for  i  =  2, ... ,  n, 

7*n  =  -^{2lQ(r,A-ri)  +  QT(r,r)UA-rt}.  (47) 

Given  a  choice  for  the  matrix  II,  Eq.  (47)  therefore  identifies  all  the  trilinear  functions  Q  for  which 
condition  (S2)  holds.  Note  that  no  restriction  is  placed  on  any  other  structural  coefficients  jijk  in 
(43)  besides  7;n ,  i  =  2, . . .  ,n.  (However,  symmetry  requires  that  71,1  =  711;  =  Jin  for  each  i.) 

Next,  consider  condition  (S3).  As  shown  by  the  foregoing  calculation,  there  exist  trilinear 
functions  Q  satisfying  (S2).  It  turns  out  the  value  of  the  expression  G(r,r,Q(r,r))  appearing  in 

(53)  does  not  depend  on  the  specific  trilinear  form  Q  we  might  choose.  This  value  is  obtained  next, 
using  conditions  (SI)  and  (S2)  as  well  as  the  results  of  Appendix  A. 

By  (SI)  and  Proposition  4,  we  have  that  Q(r,  r)  £  Es .  Thus  we  are  free  to  substitute  Q(r,  r) 
for  Aw  in  condition  (S2),  upon  which  we  directly  obtain 

3 G(r,r,Q(r,r))  =  -4lQ(r,  A~ Q(r,r))  -  2Qr (r,r)UA~ Q(r,  r) 

=  -AlQ(r,A-Q(r,r))-QT(r,r)((A-)TU  +  UA-)Q(r,r).  (48) 

Here,  A~  is  as  defined  in  Eq.  (45).  Substitution  of  (48)  into  (S3)  yields  the  following  condition 
equivalent  to  (S3): 

2 lC{r,r,r)  —  4 lQ{r,  A~ Q{r,r))  -  QT(r,r)((A~)T n  +  n A~)Q{r,r)  <  0.  (49) 

Note  that  (S3)  therefore  places  a  condition  only  on  the  quadratic  part  of  the  Liapunov  function 
candidate  V,  as  reflected  by  the  appearance  in  (49)  of  the  matrix  n. 

Conditions  (SI)  and  (S3)  are  akin  to  conditions  that  arise  in  the  stability  analysis  of  stationary 
bifurcation  for  parametrized  embeddings  of  Eq.  (1).  Under  hypothesis  (S),  such  parametrized 
systems  will  generically  exhibit  a  bifurcation  in  which  a  new  equilibrium  xe  coexists  with  the  origin 
for  each  small  |e|.  Here,  e  is  a  (normalized)  real  amplitude  parameter.  The  eigenvalue  near  zero  of 
the  bifurcated  equilibrium  xt  is  given  by  an  expansion 


20 


(3(e)  =  foe  +  (32e2  +  foe3  H - 


(50) 


To  guarantee  asymptotic  stability  of  the  new  equilibrium  xe,  one  requires  fo  =  0  and  fo  <  0.  In 
this  context  (cf.  [2]),  we  have  the  “bifurcation  formulae” 


fo  -lQ(r,r),  (51) 

fo  =  2/{C(r,  r,  r)  -  2 Q(r,  A~Q(r,  r))j.  (52) 

Note  that  condition  (SI)  is  therefore  identical  to  the  bifurcation  stability  condition  fo  =  0. 
Similarly,  (S3)  (equivalently,  Eq.  (49))  is  readily  expressed  in  terms  of  the  coefficient  fo.  Denote 
by  As(ri)  the  11-dependent  scalar 

As(ri)  :=  -QT(r,  r)((A~)T IT  +  II A~)Q(r,  r).  (53) 

Eq.  (49)  is  now  rewritten  in  terms  of  fo  and  As: 

fo  +  As(II)  <  0.  (54) 

The  following  result  is  well  known  (see,  e.g.,  [2],  [9]). 

Theorem  1.  Let  hypothesis  (S)  hold,  and  suppose  that  fo  =  0  and  fo  <  0,  where  fo  and  fo  are 
given  by  (51)  and  (52),  respectively.  Then  the  origin  of  Eq.  (1)  is  locally  asymptotically  stable. 

In  our  pursuit  of  Liapunov  functions  for  (1),  we  have  in  fact  rederived  this  result.  Indeed,  note 
that  our  condition  (SI)  requires  that  fo  —  0,  and,  for  a  given  system  (1)  for  which  fo  <  0,  the 
matrix  II  can  be  chosen  so  as  to  ensure  that  (54)  holds.  (Recall  that  II  is  any  real  symmetric  matrix 
satisfying  conditions  (i)-(iii)  of  Proposition  5.  These  conditions  are  linear  in  II.) 

An  observation  of  relevance  here  is  that,  for  any  choice  of  II,  the  quantity  A5(II)  is  nonnegative: 

AS(II)  >  0.  (55) 

Thus,  A5(n)  contributes  adversely  to  satisfaction  of  Eq.  (54). 
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We  can  now  present  the  main  result  of  this  section,  the  construction  of  a  family  of  Liapunov 

functions  V(:r)  of  the  form  (14)  for  Case  (S)  (one  zero  eigenvalue). 

Theorem  2.  Let  hypothesis  (S)  hold,  and  suppose  that  ($i  =  0  and  /?2  <  0,  where  /?i  and  /?2  are 

given  by  (51)  and  (52),  respectively.  Then  any  function  V(x)  resulting  from  Algorithm  Vs  below 

is  a  Liapunov  function  for  the  equilibrium  point  0  of  Eq.  (1). 

Algorithm  Vs.  (Construction  of  Liapunov  functions  V(ar)  =  xTVx  +  Q(x,  x,  x)  in  Case  (S)) 

Step  1.  Compute  l  and  r.  Choose  a  basis  {r2,...,rn}  for  Es .  Compute  the  dual  basis 

{l1 ,  l2, . .  . ,  /”}  to  the  basis  {r1 ,  r2, . . . ,  rn}  for  IRn .  Here,  r1  :=  r  and  /*  —  /.  Compute 

the  coefficients  f3 1  and  /?2  according  to  Eqs.  (51)  and  (52),  respectively.  Check  that  /3\  —  0 
and  /?2  <  0. 

Step  2.  Choose  any  real  symmetric  n  x  n  matrix  n  satisfying,  for  all  w  E  Es ,  w  ^  0:  (i)  nr  =  0, 

(ii)  u;rnu;  >  0,  and  (iii)  wT (A1  n  +  n A)w  <  0,  and  for  which 

|As(ri)|  <  \(32\,  (56) 

where  As(n)  is  as  defined  in  Eq.  (53). 

Step  3.  Set  V  =  l1 1  +  n. 

Step  4.  Set  the  structural  coefficients  jm,  i  =  2, . . .  ,n  to 

Tin  -  -|{2 lQ{r,  A~ A)  +  Qr(r,  r)nA-r’}.  (57) 

Step  5.  Symmetry  requires  that  7^  be  independent  of  permutations  in  the  indices  i,j,k.  The 
structural  coefficients  in  the  representation 

g(x\x2,x3)  =  EEE  7ijk{lix1)(Px*)(lkx3)  (58) 

i=l j= 1 k= 1 

which  have  not  been  specfied  in  Steps  1-4  are  either  determined  by  symmetry  and  Eq.  (57), 
or  can  be  assigned  arbitrarily,  subject  only  to  the  symmetry  requirement. 
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7  CONSTRUCTION  OF  LIAPUNOV  FUNCTIONS  IN  THE 


CASE  OF  A  PAIR  OF  PURE  IMAGINARY  EIGENVALUES 

In  this  section,  Liapunov  functions  are  constructed  for  the  origin  of  Eq.  (1)  under  hypothesis  (H). 
The  construction  parallels  that  of  the  last  section,  while  differing  from  it  in  several  respects.  For 
example,  the  Liapunov  function  candidates  used  in  the  previous  section  consist  only  of  quadratic 
and  cubic  terms  in  the  state.  In  this  section,  quartic  terms  also  appear  in  the  Liapunov  function 
candidates  (cf.  Eq.  (15)).  If  quartic  terms  were  not  included  in  the  assumed  form  of  the  Liapunov 
function  candidates,  the  construction  would  fail  generically  in  Case  (H).  A  second  difference  be¬ 
tween  the  calculations  of  this  and  the  preceding  section  concerns  the  adoption  of  complex  notation. 
Although  not  essential,  complex  notation  is  both  natural  and  convenient  when  considering  local 
asymptotic  stability  of  systems  (1)  under  hypothesis  (II). 

7.1  Conditions  for  V  <  0 

Consider,  then,  the  time-derivative  of  the  Liapunov  function  candidate 

V(x)  =  xTVx  +  Q{x,  x,  x)  +  T(x,  x,  x,  x)  (59) 

along  trajectories  of  Eq.  (1).  Using  Proposition  1,  and  the  Taylor  series  representation  (13)  of 
/(x),  we  find  that  this  derivative  is  given  by 

V(x)  =  xr  (A1  V  +  VA)x 

+  2QT(x,  x)Vx  +  3 <?(x,  x,  Ax) 

+  2 CT(x,  x,x)Vx  +  3f/(x,  x,  Q(x,  x))  +  4T (x,  x,  x,  Ax)  +  •  •  •  (60) 

In  Section  5,  it  was  noted  that  any  vector  x  E  IRn  has  a  unique  representation 

x  =  ar  +  ar  +  w.  (61) 

Here,  r  is  the  right  eigenvector  of  A  associated  with  the  eigenvalue  iojc  (as  specified  in  Section  1),  f 
is  the  complex  conjugate  of  r  and  as  such  is  a  right  eigenvector  of  A  associated  with  the  eigenvalue 
— iu>c ,  a  is  a  complex  scalar,  and  w  E  Es . 

The  first  step  in  our  procedure  for  obtaining  conditions  for  local  negative  definiteness  of  V(x) 
is  to  substitute  in  Eq.  (60)  the  representation  (61)  for  x  and  the  representation  (37)  for  the  matrix 
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V,  and  group  terms  according  to  their  order  in  |(a,a,u;)|.  (Recall  that  (37)  states  that,  in  Case 
H,  V  is  chosen  from  among  matrices  of  the  form  V  =  1 T  +  lTI+ l  Tl,  where  II  satisfies  conditions 
(i)-(iii)  of  Proposition  5.)  Using  Proposition  4  and  Remark  2  (orthogonality  of  left  and  right 
eigenvectors),  and  after  a  considerable  amount  of  algebra  and  reordering  of  terms,  we  can  obtain 
explicit  formulae  for  the  quadratic,  cubic  and  quartic  terms  in  |(a,a,tn)|  (these  are  [V](2p  [V](3) , 
and  [V](4),  respectively)  in  the  expansion 

V(x)  :=  [V](2)  +  [V](3)  +  [V](4)  H -  (62) 

of  V(*).  (This  is  in  analogy  with  Eqs.  (39)-(42)  of  the  preceding  section.)  For  example,  one  can 
check  that  [V](2j,  [V](3\  are  given  by 

=  wt{At  n  +  UA)w,  (63) 

=  ci3{2lQ(r,r)  +  3iu;cG(r,r,r)} 

+  a3 {2 lQ(r,  r)  -  3 iu>cG(r,  r,  r)} 

+  a2a{2/Q(r,  r)  +  4 lQ(r,  r )  +  3 iu>c<7(r,  r,  f)} 

+  a2a{2lQ(r,  f)  +  4 lQ(r,  r)  -  3 iu)cQ{r,  r,  r)} 

+  a2{2 Q7  (r,  r)IIu;  +  4 lQ(r,  w)  +  3 Q{r,  r,  Aw)  +  6iu>cQ(r ,  r,  in)} 

+  a2{2QT(r,  r)IIu;  +  4/Q(f ,  in)  +  3 (7(f,  f,  Aw)  —  6 iuicQ(r,  f,  in)} 

+  aa{4QT(r,  f)IItn  +  4/Q(r,  w)  +  4 lQ(r,  in)  +  QQ(r,  r,  Ain)} 

+  a{2/Q(in,  in)  +  6 Q(r,  in,  Ain)  +  3 iwcQ{r,  in,  in)} 

+  a{2/Q(in,  in)  +  6 Q(r,  in,  Ain)  —  3iwcC/(f,  in,  in)} 

+  ZQ[w,w,  Ain).  (64) 

Note  that  both  [V]^)  and  [V]^  as  given  above  are  real-valued,  as  expected.  There  is  no  need  to 
give  the  full  expression  for  [V]^  here.  Instead,  we  proceed  directly  to  the  statement  of  the  following 
preliminary  result  in  the  construction  of  Liapunov  functions  under  hypothesis  (II).  In  the  proof  of 
this  result,  which  is  analogous  to  Proposition  6  in  the  preceding  section,  values  of  certain  pertinent 
terms  appearing  in  the  expansion  of  [V](4j  will  be  given.  All  the  terms  in  the  expansion  of  [V]^ 
may  be  obtained  readily,  to  result  in  an  expression  analogous  to  (64). 

Proposition  7.  Let  hypothesis  (H)  hold.  Suppose  that,  for  some  real  matrix  II  satisfying  the 


[V](2) 

[V](3) 
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conditions  of  Proposition  5,  symmetric  real  trilinear  function  Q(x1 ,  x2 ,  z3),  and  symmetric  real 
tetralinear  function  T(ad ,  x2 1  x3,  a;4),  the  following  seven  conditions  hold: 

(HI)  2 [Q(r,  r)  +  3iwc(7(r,  r,  r)  =  0, 

(H2)  2 lQ(r,  r )  +  4 lQ(r,  r)  +  3 iu>cQ(r,  r,  r )  =  0, 

(H3)  2QT(r,  r)IIu)  +  4/Q(r,  ui)  +  3 Q(r,  r,  ( A  +  2iwcI)w)  —  0,  for  all  w  £  Es, 

(H4)  2 Qr(r,  r)IIu;  +  2/Q(r,  w )  +  2/Q(r,  w)  +  3 G(r,  r,  Aw)  —  0,  for  all  w  £  Es , 

(H5)  2lC(r,  r,  r)  +  3 Q(r,  r,  Q(r,  r ))  +  4iwcT (r,  r,  r,  r)  =  0, 

(H6)  IC(r,  r,  r)  +  3/C(r,  r,  r)  +  +3^(r,  r,  Q(r,  f))  +  3 Q(r,  r,  Q(r,  r))  +  4iwcT (r,  r,  r,  r)  =  0, 

(H7)  3£e{2 /C(r,  r,  r)  +  t/(r,  r,  Q(r,  f))}  +  2 £/(r,  r,  Q(r,  r))  <  0. 

Then  the  origin  of  Eq.  (1)  is  locally  asymptotically  stable. 

Proof.  The  proof  consists  of  a  judicious  application  of  Lemma  1,  which  gives  general  sufficient 
conditions  for  local  negative  definiteness  of  a  class  of  bivariate  functions.  Identify  the  scalar  variables 
u  and  v  in  Lemma  1  as  u  :=  |«j|,  v  :=  |a|  =  |a|.  The  bivariate  function  S(u,  v)  of  Lemma  1  is  taken 
to  be  a  local  upper  bound  on  V(a:).  We  show  that  conditions  (H1)-(H7)  are  sufficient  for  there 
to  exist  such  an  upper  bound  6(u,v)  satisfying  the  hypotheses  of  Lemma  1.  More  specifically,  we 
show  that  (H1)-(H7)  results  in  satisfaction  of  the  simplified  requirements  of  Remark  1(a)  following 
Lemma  1  (this  adds  the  condition  ai2  =  0).  To  ensure  absence  of  the  u3-term  in  6(u,v),  i.e.,  that 
ao3  =  0,  we  require  the  coefficients  of  a3,  a3,  a2a  and  a2a  in  Eq.  (64)  to  vanish.  The  coefficient  of 
a3  is  precisely  the  expression  on  the  left  side  of  (HI).  (Note  that  the  a3-term  in  (64)  is  the  complex 
conjugate  of  the  a3-term,  and  thus  also  vanishes  when  (HI)  is  in  force;  analogous  comments  apply 
below.)  Similarly,  the  expression  on  the  left  side  of  (H2)  is  simply  the  coefficient  of  a2d  in  Eq. 
(64).  Thus,  (HI)  and  (H2)  combined  ensure  that  o03  =  0.  To  ensure  absence  of  the  ur2-term  in 
8(u,v),  i.e.,  that  ai2  =  0,  we  require  the  “linear-in-in”  coefficients  of  a2  and  ad  in  Eq.  (64)  to 
vanish  for  each  w  £  Es .  (The  coefficient  of  a2,  being  the  conjugate  of  that  of  a2,  will  then  vanish 
automatically.)  Inspection  of  Eq.  (64)  reveals  that  this  is  equivalent  to  conditions  (H3)  and  (114) 
above.  It  remains  to  show  that  020  <  0  and  004  <  0  (in  the  notation  of  Lemma  1).  That  020  <  0 
follows  immediately  from  Eq.  (63)  and  Corollary  1.  Conditions  ensuring  that  a0 4  <  0  can  only 
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result  from  examination  of  the  quartic  terms  in  V.  However,  quartic  terms  in  [V]^4^  which  involve  w 
are  irrelevant  to  this  requirement.  The  needed  coefficients  can  be  obtained  readily  by  substituting 
Eq.  (61)  for  x  in  the  formula 

[V](4)(x)  =  2 CT(x,  x,  x)Vx  +  3 Q(x,x,Q(x,  x))  +  4 T(x,x,  x,Ax),  (65) 

using  the  orthogonality  of  left  and  right  eigenvectors  (cf.  Proposition  4  and  Remark  2),  and  using 
the  fact  that  a  multilinear  form  is  linear  in  each  argument.  The  expansion  of  [V](4)(ur  -f  dr  +  w ) 
is  seen  to  contain  five  terms  which  are  quartic  in  (a,  a):  an  a4-term,  an  a3a-term,  the  conjugates 
of  these  two,  and  an  a2a2-term.  Of  these  five,  only  the  latter  is  sign-definite:  a2a2  =  |a|4.  Thus, 
we  require  the  coefficient  of  a2  a2  in  the  expansion  of  [V]j4j(ar  +  ar  +  w)  to  be  negative,  and  the 
coefficients  of  a4  and  a3a  to  vanish.  It  is  readily  verified  that  the  left  side  of  (H5)  is  the  coefficient 
of  a4,  and  the  left  side  of  (H6)  is  half  the  coefficient  of  a3d.  Finally,  it  is  straightforward  to  check 
that  the  left  side  of  (H7)  is  one-sixth  the  coefficient  of  a2a 2  in  the  expansion  of  [V]^4j(ar  +  dr  +  u>). 
Thus,  conditions  (H5)-(H7)  together  give  the  desired  negativity  of  ao4-  ■ 

7.2  Algorithm  for  Construction  of  V  in  Case  (H) 

Conditions  (H1)-(H7)  may  be  solved  for  a  trilinear  function  f7(x4,  x2 ,  x3)  and  a  tetralinear  function 
T(x4,  x2,  x3,  x4),  under  an  appropriate  auxiliary  condition  ensuring  stability  of  the  origin  and  for 
a  given  matrix  H  satisfying  the  conditions  of  Proposition  5.  The  procedure  is  much  the  same  as 
was  carried  out  in  the  preceding  section,  where  the  coordinate  representation  of  a  trilinear  form 
£7(x4,x2,x3)  was  employed  to  solve  (S1)-(S3)  for  the  structural  coefficients  jtJi.  of  Q.  Due  to  this 
similarity,  only  a  summary  of  the  main  steps  in  the  derivation  is  deemed  necessary  here,  with  the 
result  for  the  Liapunov  functions  we  obtain  summarized  below  in  Algorithm  Vh- 

For  convenience,  we  continue  to  employ  complex  notation,  and  choose  a  coordinate  basis 
{r1 ,  r2, . . . ,  r”}  for  IRn  in  which  r1  :=  r,  r2  :=  r,  and  r3, . . . ,  rn  lie  in  Es  C  IRn ■  To  this  ba¬ 
sis  there  corresponds  a  unique  dual  basis  of  row  vectors  j/1 ,  l2, . . . ,  ln},  where  l 1  :=  l  and  l'2  :=  l. 
The  trilinear  and  tetralinear  functions  Q^x1 ,  x2 ,  x3)  and  T(x4,  x2,  x3,  x4)  are  then  expressed  in  the 
coordinate  representations  (17)  and  (18),  respectively.  We  seek  the  minimum  set  of  specifications  on 
the  associated  structural  coefficients  jijk  and  T{jkp,  respectively,  under  which  conditions  (H1)-(II7) 
above  hold. 
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Conditions  (HI)  and  (H2)  are  interpreted  in  this  framework  simply  as  assigning  values  to  the 
structural  coefficients  7ni  and  7112,  respectively.  Next  consider  (H3).  Since  Es  is  invariant  under 
A,  w  £  Es  implies  that  ( A-\-2iucI)w  lies  in  the  complexification  of  Es  (also  referred  to  as  Es  below). 
Thus,  we  can  define  a  vector  w  :=  ( A  +  2iucI)w,  noting  that  the  matrix  inverse  in  the  equation 
w  =  (A  4-  2iu>cI)~1w  exists  by  hypothesis  (H).  Interpreting  (H3)  as  a  requirement  on  each  basis 
vector  r3, . . .  ,rn  of  Es ,  we  find  that  (H3)  amounts  to  a  specification  of  the  structural  coefficients 
jin,  i  —  3,4  ,...,n.  Similarly,  (H4)  amounts  to  a  specification  of  the  structural  coefficients  Jn2, 
i  —  3, 4, ... ,  n. 

Since  each  of  the  structural  coefficients  7,11,  i  =  l,2,...,n  is  fixed  (by  one  of  (II1)-(H3)), 
Q(r ,  r,  x )  is  fixed  for  any  x  £  <En .  The  coefficients  JH2,  i  =  1,2, ...  ,n  are  also  fixed:  7^2  is  fixed  by 
(H2),  7,12,  *  =  3,4,.. .  ,n  are  fixed  by  (114),  and  we  also  have  7212  =  7112  by  applying  Proposition 
3  to  (H2).  Thus,  Q(r,r,x )  is  also  fixed  for  any  x  £  <Vn .  By  these  remarks,  it  follows  that  the 
terms  G(r,r,Q(r,r)),  Q(r,  r,  Q(r,  fi)),  G(r,  r,  Q(r,  r))  appearing  in  (H5)  and  (H6)  are  determined  by 
(H1)-(II4).  Their  values  may  be  found  by  expressing  Q(r,  r)  and  Q(r,r )  as  linear  combinations  of 
the  basis  vectors  r*,  i  =  1  ,...,n,  and  then  employing  the  coordinate  representation  of  Q.  Thus, 
(H5)  and  (H6)  serve  to  assign  the  values  of  T(r,r,r,r)  (=  run)  and  T(r,r,r,r)  (=  77112)-  The 
importance  of  including  the  quartic  term  T  in  the  Liapunov  function  candidate  now  becomes  clear: 
with  T  =  0,  (H5)  and  (H6)  become  constraints  on  the  system  which  do  not  constitute  necessary 
conditions  for  stability.  However,  with  inclusion  of  a  quartic  term  T,  (H5)  and  (H6)  are  quite  easily 
satisfied. 

By  the  remarks  above,  it  follows  that  the  quantity  appearing  on  the  left  side  of  (117)  is  completely 
specified  by  the  system  dynamics  and  the  matrix  n.  Next,  we  sketch  the  derivation  of  an  explicit 
reformulation  of  the  left  side  of  (H7)  in  terms  of  system  (1)  and  n.  A  stability  coefficient  which 
arises  in  the  study  of  Hopf  bifurcation  for  parametrized  versions  of  (1)  under  hypothesis  (IT)  will 
appear  in  the  reformulation.  The  value  of  this  coefficient,  which  we  denote  as  /?2,  is  recalled  next. 
Note  that  this  coefficient  /?2  is  distinct  from  the  coefficient  of  the  same  name  appearing  in  Section 
6.  In  the  present  context,  /?2  relates  to  an  (even  in  e)  expansion  /?(e)  =  ^e2  +  /^e4  +  •  •  •  of  the 
Floquet  exponent  near  zero  of  bifurcated  periodic  solutions  of  parametrized  embeddings  of  Eq.  (1). 
(Compare  with  Eq.  (50)  for  the  analogous  eigenvalue  expansion  in  Case  (S).) 
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Define  vectors  £  and  77  by 


£  :=  “2^  1Q(rd),  (66) 

77  :=  I(2io;c/-^)-:1g(r,r).  (67) 

(In  [9]  and  [1],  £  and  77  are  denoted  as  a  and  6,  respectively.)  Then  /?2  is  given  by  the  “bifurcation 
formula”  [9],  [1] 

P 2  :=  2*e{2/Q(r,0  +  /Q(r,  77)  +  |/C(r,  r,  f)}.  (68) 

Conditions  (H1)-(H4)  can  be  used  to  replace  (H7)  with  an  equivalent  condition  stated  explicitly 
in  terms  of  the  Taylor  expansion  (13)  of  system  (1).  Consider  the  two  terms  G{r,r,Q(r,r))  and 
G(r,  r,  Q(r,  r))  appearing  in  (H7).  The  former  quantity  is  of  the  form  Q(r,  r,  x),  which  occurs  in  the 
statements  of  conditions  (HI)  (with  x  =  r),  (H2)  (with  x  =  r)  and  (113)  (with  x  G  Ea).  Thus,  we 
resolve  the  vector  x  (=  Q(r ,  r))  into  its  components  lying  in  the  subspace  Es  and  along  the  r-  and 
f-directions,  and  then  apply  (H1)-(H3).  The  r-component  of  any  vector  x  is  given  by  ( lx)r  (recall 
the  normalization  lr  =  1),  the  r-component  is  (lx)r,  and  the  ^-component  is  the  remainder 


xs  :=  x  —  ( lx)r  —  ( lx)r .  (69) 

Using  (II1)-(H3),  we  have  that,  for  any  x, 

G(r,r,x)  -  -^—(lx)lQ(r,r)-^—(lx)(lQ(r,r)  +  2jQ(r,r)) 

-  ^{2 QT(r,  r)n(^l  +  2iwc/)"V  +  4 lQ(r,  (A  +  2 iucI)~lxs}.  (70) 

Since  r  is  an  eigenvector  of  A,  it  is  also  an  eigenvector  of  ( A  +  2iu>cI)~1  and  of  {A  —  2  iuj)-1. 
Using  this  fact,  and  the  fact  that  nr  =  0,  this  expression  may  be  expanded  and  simplified.  Letting 
x  =  Q(f,f),  the  resulting  expression  is 

U(r,r,  Q(r,f))  =  ^/Q(r,  77)  +  ^QT(r,  r)n?7 

“  r))(lQ(r,  r))  -  ^-(/Q(r,  r))(/Q(r,  r)).  (71) 
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Thus 


5fte(£(r,  r,  x))  =  $te{^IQ(r,  r?)  +  ^QT(r,  r)II^}.  (72) 

The  following  formula  for  <?(r,  r,  <2(r,  r))  is  obtained  in  a  similar  fashion.  The  fact  that  r  and  f 
are  eigenvectors  of  A~l  is  employed  in  the  computation. 


G(r,  f,  Q(r,  f))  =  |»c{/Q(r,0  +  ^QT(r,  r)n^}.  (73) 

Condition  (H7)  may  now  be  rewritten  explicitly  as 


ft  +  Atf(II)  <  0,  (74) 

where  /?2  is  as  defined  above,  and  where  A h  is  given  by  the  real  number 

Atf(II)  :=  Q^ffdA-'fn  +  nA-^Q^r) 

~  \qT{ r,  r){n(T  +  2icuc7)-1  +  ((A  -  2iu;ciyl)TU}Q(f)  r).  (75) 

It  is  not  difficult  to  ascertain  that  A h  is  nonnegative,  although  this  is  not  an  essential  consider¬ 
ation.  (One  proof  of  this  uses  the  fact  that  formula  (75)  is  a  special  case  of  the  sum  of  a  Hermitian 
form  and  a  quadratic  form.) 

We  have  just  rederived  the  following  known  criterion  for  asymptotic  stability  in  Case  (II). 

Theorem  3.  Let  hypothesis  (II)  hold,  and  suppose  that  /?2  <  0,  where  /?2  is  given  by  Eq.  (68). 
Then  the  origin  of  (1)  is  locally  asymptotically  stable. 

Regarding  the  structural  coefficients  ~fijk  and  T,-jfcp  that  have  not  been  specified  explicitly  in  the 
foregoing  analysis,  only  two  constraints  remain:  The  first,  the  symmetry  requirement,  entails  that 
the  value  of  a  coefficient  is  independent  of  the  order  of  subscript  indices.  The  second  constraint 
is  that  the  function  V(x)  must  be  real-valued  for  x  £  IRn .  Using  Proposition  3,  it  is  found  that 
this  latter  requirement  is  equivalent  to  what  might  be  called  a  conjugate  symmetry  relationship 
among  the  structural  coefficients,  the  exact  nature  of  which  is  specified  in  the  next  Corollary  to 
Proposition  3. 
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Corollary  2.  (Conjugate  Symmetry  of  Structural  Coefficients  in  Case  (H))  Denote,  for 
any  positive  integer  i,  the  quantity  [i]  (the  “complement  of  i”) 


W  := 


2  if  *  =  1, 

1  if  i  —  2, 


(76) 


l  i  otherwise 

Then,  V(:r)  is  real- valued  for  each  x  £  IRn  if  and  only  if  the  structural  coefficients  7,-yj.  and  77^ 
satisfy  the  following  relationship: 


7ijk  ~  7[t][,p]>  (77) 

Tijkp  =  r[i][7][fc][p]-  (78) 

Proof.  Follows  immediately  from  Proposition  3.  ■ 

The  foregoing  construction  of  a  family  of  Liapunov  functions  V(x)  of  the  form  (59)  for  the  case 

in  which  Df( 0)  possesses  a  pair  of  pure  imaginary  eigenvalues  is  summarized  in  the  next  result  and 
algorithm. 


Theorem  4.  Let  hypothesis  (H)  hold,  and  suppose  that  <  0,  where  /%  is  given  by  Eq.  (68). 

Then  any  function  V(:r)  resulting  from  Algorithm  Vh  below  is  a  Liapunov  function  for  the  equilib¬ 
rium  point  0  of  Eq.  (1). 

Algorithm  V//.  (Construction  of  Liapunov  functions  V(:c)  =  xTVx  +  Q(x,  x,  x)  +  T(x,  x,  x,  x)  in 

Case  (H)) 

Step  1.  Compute  /  and  r  and  choose  a  basis  {r1,  r2, . . . ,  r”)  for  IRn  with  r1  :=  r,  r2  :=  r,  and  for 
which  {r  , . . .  ,r")  is  a  basis  for  Es .  Calculate  the  row  vectors  {/3, . . . ,  ln}  of  the  associated 
dual  basis.  Check  that  /%  <  0,  where  /?2  is  as  defined  in  (68). 

Step  2.  Pick  any  II  satisfying  (i)  Ilr  =  Ilr  =  0,  (ii)  ru^IIu;  >  0,  and  (iii)  wT(ATH  +  UA)w  <  0  for 
all  w  £  Es,w  ^  0,  and  such  that  A#(II)  <  |/?2 1 ,  where  A//(II)  is  given  by  Eq.  (75). 

Step  3.  Set  V  —  II  +  lTl  -(- 1  Tl. 

Step  4.  Set  Q(x1 ,  x2,  x3)  =  Ef=i  Ej=i  Efc=i  HjkQ'x1)^3  x2)(lkx3),  and 

Tin  =  “3 £“^Kr>r)>  (79) 
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7112 


=  -^m(r,r)  +  2lQ(r,r)},  (80) 

Till  =  ~{4lQ{r,(2iu>cI  +  A)-\i)  +  2QT(r,r)Il(A  +  2iucI)-1ri},  i  =  3, . . .  ,  n,  (81) 
7*12  =  -^{TQ(r,A~1rt)  +  lQ(r,A~1r,)  +  QT(r,f)UA~1r1},  i  =  3, . . . ,  n.  (82) 

Step  5.  Set  T(xl,x2,  x3,  x4)  =  E"=i  E"=  i  Efc=i  Ep=i  Tijkp{llxx){P x2){lkxz){PxA).  Here,  rmi  (= 
T(r,  r,  r,  r))  and  rm2  (=  T(r,  r,  r,  r))  are  selected  according  to  (H5)  and  (H6),  respectively. 

Step  6.  All  structural  coefficients  7 ijk  and  r,jfcp  which  have  not  been  specified  in  Steps  1-5  are 
assigned  arbitrarily,  modulo  the  symmetry  requirement  and  the  conjugate  symmetry  require¬ 
ment  (Eqs.  (77),  (78)). 


8  CONCLUDING  REMARKS 

Liapunov  functions  for  nonlinear  systems  with  either  of  the  two  simplest  critical  cases  have  been 
explicitly  constructed.  For  the  case  in  which  the  sytem  linearization  possesses  a  simple  zero  eigen¬ 
value,  generically  the  Liapunov  functions  need  contain  only  quadratic  and  cubic  terms  in  the  state. 
However,  when  a  complex  conjugate  pair  of  simple,  pure  imaginary  eigenvalues  are  present,  the 
Liapunov  functions  contain  quartic  terms,  in  addition  to  the  quadratic  and  cubic  terms.  These 
Liapunov  functions  were  shown  to  predict  local  asymptotic  stability  precisely  when  certain  known 
sufficient  conditions  from  bifurcation  analysis  are  satisfied.  We  have  obtained  parametrized  “fam¬ 
ilies”  of  Liapunov  functions  for  the  studied  critical  cases,  in  the  same  sense  that  the  Liapunov 
matrix  equation  yields  an  infinite  set  of  quadratic  Liapunov  functions  for  asymptotically  stable  lin¬ 
ear  time-invariant  systems.  The  Liapunov  functions  are  computed  directly  in  terms  of  the  Taylor 
series  expansion  of  the  vector  field  f(x),  and  are  thus  amenable  to  symbolic  computer  coding.  The 
use  of  these  Liapunov  functions  in  the  design  of  feedback  control  laws  for  critical  nonlinear  systems 
is  a  topic  for  future  investigation. 
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APPENDIX  A.  SOLUTION  OF  LINEAR  ALGEBRAIC  EQUATIONS 
WITH  SINGULAR  COEFFICIENT  MATRIX 

Consider  the  system  of  linear  equations 

Ax  —  b  (A.l) 

where  A  is  a  real  n  x  n  matrix  and  b  £  IRn.  Suppose  that  A  has  a  simple  zero  eigenvalue.  Let 
r  and  /  denote  right  (column)  and  left  (row)  eigenvectors  of  A,  respectively,  corresponding  to  the 
zero  eigenvalue,  and  require  that  these  be  chosen  to  satisfy  Ir  =  1.  Under  these  conditions,  the 
Fredholm  Alternative  asserts  that  (A.l)  has  a  solution  if  and  only  if  lb  =  0.  Moreover,  the  Fredholm 
Alternative  also  implies  that,  if  (A.l)  has  a  solution  x°,  then  the  totality  of  solutions  is  given  by 
the  one-parameter  family  x  =  x°  +  ar  where  a  £  IR  is  arbitrary.  The  solution  is  rendered  unique 
upon  imposing  a  normalization  condition  which  specifies  the  value  of  lx. 

Introduce  subspaces  Ec,  Es  C  IRn  as  follows:  Ec  is  the  one-dimensional  subspace 

Ec  :=  span{r},  (A. 2) 

and  Es  is  the  (n  —  l)-dimensional  subspace 

Es  :=  {x  £  IRn\  lx  =  0}.  (A. 3) 

From  the  foregoing,  we  have  in  particular  that  if  lb  =  0  then  the  system  Ax  =  b,  lx  =  0  has  a 
unique  solution.  Equivalently,  (A.l)  has  a  unique  solution  in  Es  for  any  vector  b  £  Es .  This  proves 
that  the  restriction  [8,  p.  199]  A \E,  of  the  linear  map  A  to  Es  defines  an  invertible  (one-to-one  and 
onto)  map.  In  the  next  result,  we  exhibit  the  unique  solution  which  lies  in  Esof  the  system  Ax  =  b, 
lx  =  0.  The  proof  is  elementary  [2]. 

Proposition  A.l  The  unique  solution  of  Ax  —  b,  lx  =  0  given  that  lb  —  0  is 

x-(ATA  +  lTl)~1ATb.  (A. 4) 

This  result  motivates  the  following  introduction  of  notation: 

A~  :=  (ATA  +  lTl)~lAT.  (A. 5) 

Thus,  the  inverse  of  the  restricted  map  A\Es  exists  and  is  given  by 

(A\es)~'=A-.  (A. 6) 


32 


ACKNOWLEDGMENT 

The  authors  are  grateful  to  Dr.  Der-Cherng  Liaw  for  helpful  comments.  They  also  thank  Reviewers 
3  and  4  for  their  suggestions.  This  research  has  been  supported  in  part  by  the  National  Science 
Foundation’s  Engineering  Research  Centers  Program:  NSFD  CDR-88-03012,  by  the  Air  Force  Office 
of  Scientific  Research  under  URI  Grant  AFOSR-90-0015,  by  NSF  Grant  ECS-86-57561,  and  by  the 
TRW  Foundation. 

REFERENCES 

[1]  E.H.  Abed  and  J.-H.  Fu,  “Local  feedback  stabilization  and  bifurcation  control,  I.  Hopf  bifur¬ 

cation,”  Systems  and  Control  Letters  Vol.  7  (1986),  pp.  11-17. 

[2]  E.H.  Abed  and  J.-H.  Fu,  “Local  feedback  stabilization  and  bifurcation  control,  II.  Stationary 

bifurcation,”  Systems  and  Control  Letters  Vol.  8  (1987),  pp.  467-473. 

[3]  D.J.  Allwright  and  A.I.  Mees,  “Lyapunov  method  applied  to  the  Hopf  bifurcation,”  J.  London 

Math.  Soc.  (2)  20  (1979),  pp.  293-299. 

[4]  R.L.  Bishop  and  S.I.  Goldberg,  Tensor  Analysis  on  Manifolds,  Dover  Publications,  New  York, 

1980. 

[5]  J.  Carr,  Applications  of  Centre  Manifold  Theory,  Springer-Verlag,  New  York,  1981. 

[6]  S.N.  Chow  and  J.K.  Hale,  Methods  of  Bifurcation  Theory,  Springer-Verlag,  New  York,  1982. 

[7]  J.-H.  Fu,  “Liapunov  functions  and  stability  criteria  for  nonlinear  systems  with  multiple  critical 

eigenvalues,”  manuscript  in  preparation. 

[8]  K.  Hoffman  and  R.  Kunze,  Linear  Algebra,  Second  Edition,  Prentice-IIall,  Englewood  Cliffs, 

N.J.,  1971. 

[9]  L.N.  Howard,  “Nonlinear  oscillations,”  in:  F.C.  Hoppensteadt,  Ed.,  Nonlinear  Oscillations  in 

Biology  (Amer.  Math.  Society,  Providence,  RI,  1979)  pp.  1-67. 

[10]  G.  Iooss  and  D.D.  Joseph,  Elementary  Stability  and  Bifurcation  Theory,  Springer-Verlag, 
New  York,  1976. 


33 


[11]  N.N.  Krasovskii,  Stability  of  Motion,  Stanford  Univ.  Press,  Stanford,  CA,  1963.  (Translation, 
by  J.L.  Brenner,  of  the  1959  Russian  edition.) 

[12]  H.G.  Kwatny  and  X.-M.  Yu,  “Energy  analysis  of  load-induced  flutter  instability  in  classical 
models  of  electric  power  networks,”  IEEE  Trans.  Circuits  and  Systems,  Vol.  36,  pp.  1544- 
1557,  1989. 

[13]  J.P.  LaSalle  and  S.  Lefschetz,  Stability  by  Liapunov’s  Direct  Method  with  Applications, 
Academic  Press,  New  York,  1961. 

[14]  A.M.  Lyapunov,  “Probleme  general  de  la  stabilite  du  mouvement,”  Ann.  Fac.  Sci.  Toulouse, 
Vol.  9,  pp.  203-474,  1907.  Reprinted  in  Ann.  Math.  Studies,  No.  17,  Princeton  Univ.  Press, 
Princeton,  NJ,  1947. 

[15]  I.  Malkin,  “Some  basic  theorems  of  the  stability  of  motion  in  the  critical  cases,”  AMS  Trans¬ 
lations  38,  1958. 

[16]  M.  Marcus,  Finite  Dimensional  Multilinear  Algebra,  (Part  I),  Marcel  Dekker,  New  York,  1973. 

[17]  A.I.  Mees  and  L.O.  Chua,  “The  Hopf  bifurcation  theorem  and  its  applications  to  nonlinear 
oscillations  in  circuits  and  systems,”  IEEE  Trans.  Circuits  and  Systems  CAS-26,  1979,  pp. 
235-254. 

[18]  C.E.  Pearson,  Ed.,  Handbook  of  Applied  Mathematics,  Second  Edition,  Van  Nostrand  Rein¬ 
hold,  New  York,  1983. 

[19]  E.D.  Sontag,  “Feedback  stabilization  of  nonlinear  systems,”  pp.  61-81  in:  Robust  Control 
of  Linear  Systems  and  Nonlinear  Control,  Vol.  4,  M.A.  Kaashoek,  J.H.  van  Schuppen  and 
A.C.M.  Ran,  Eds.,  Birkhauser,  Boston,  1990. 

[20]  V.I.  Zubov,  Methods  of  A.M.  Lyapunov  and  their  Application,  P.  Noordhoff,  Ltd.,  Groningen, 
1964.  (English  translation  under  the  auspices  of  the  U.S.  Atomic  Energy  Commission.  L.F. 
Boron,  Editor). 


34 


